	//#ifndef LAPLACE2_HPP
	//#define LAPLACE2_HPP
#include "Erf.hpp"
#include <armadillo>
#ifndef MAXIT
#define MAXIT 1000
#endif
using namespace arma;

#define _THRES 0.0000001

class Laplace
{
	public:

	Laplace(void)
	{	
		boost::math::normal_distribution<> *d=new boost::math::normal_distribution<>(0,1); 
		_d=d;
		_conv=0;
	}
	~Laplace()
	{
		delete _d;
	}
	void operator()(mat X,double *Y,mat m,double si,mat init)
	{
		double e=1;
		double thres=_THRES;
		int p=X.n_cols;
		mat theta(p,1);
		mat Ym=vect2mat<double>(Y,X.n_rows,1);
		if(init(0,0)==0)
		{
			theta=inv(X.t()*X)*X.t()*Ym;
		}else{
			theta=init;
		}
		mat b(p,1);
		int i=0;
		while(e>thres & MAXIT>i)
		{
			b=theta-inv(Hessian(X,Y,theta,si))*Gradiant(X,Y,theta,si);
			e=abserr(theta-b);	
			theta=b;
			i++;
		}
		if(i>MAXIT){cout << "reached max iteration!";_conv=1;	}
		_m=theta;
		_S=-inv(Hessian(X,Y,theta,si));
		

	}

	inline double abserr(mat X)
	{
		double sum=0;
		int n=X.n_rows;
		for(int i=0;i<n;i++)
		{
			sum+=abs<double>(X(i,0));
		}
		return sum;
	}

	inline mat Gradiant(mat X,double *Y,mat theta,double si)
	{
		int p=X.n_cols;
		mat sum(p,1);
		sum.fill(0);
		int n=X.n_rows;
		double xb=0;
		for(int i=0;i<n;i++)
		{
			xb=as_scalar(X(i,span::all)*theta);
			double pxb=pdf(*_d,xb);
			double gpxb=cdf(*_d,xb);
				
			if(Y[i]==1)
			{	
				sum=sum+pxb*(1/gpxb)*X(i,span::all).t();				
			}else{
							
				sum=sum-pxb*(1/(1-gpxb))*X(i,span::all).t();				
			}
		}
		return sum-as_scalar(arma::sum(theta,0))/si;
		
	}

	inline	mat Hessian(mat X,double *Y,mat theta,double si)
	{
		int p=X.n_cols;
		mat sum(p,p);
		sum.fill(0);
		int n=X.n_rows;
		double xb=0;
		for(int i=0;i<n;i++)
		{
			
			xb=as_scalar(X(i,span::all)*theta);
			double pxb=pdf(*_d,xb);
			double gpxb=cdf(*_d,xb);
			mat xx=X(i,span::all).t()*X(i,span::all);	
			if(Y[i]==1)
			{	
				sum=sum+(xb*pxb/gpxb+pow(pxb,2)/pow(gpxb,2))*xx;				
			}else{
				sum=sum+(-xb*pxb/(1-gpxb)+pow(pxb,2)/pow(1-gpxb,2))*xx;				
			}
		}
		return -sum-(double)1/si;

	}
	
	mat Get_Sig(void){return _S;}
	mat Get_Mu(void){return _m;}
	int Get_conv(void){return _conv;}
	private:
	mat _m;
	mat _S;
	int _conv;
	boost::math::normal_distribution<> *_d; 
};
#include "LALogit.hpp"
#include "VB.hpp"
#include "LaplaceLS.hpp"
//#endif
